library(tidyverse)
library(rpart)Homework 4
QTM 220
\[ \DeclareMathOperator{\E}{E} \DeclareMathOperator{\CATE}{CATE} \DeclareMathOperator{\Var}{V} \DeclareMathOperator{\hVar}{\widehat{V}} \DeclareMathOperator{\bias}{bias} \DeclareMathOperator*{\argmin}{argmin} \newcommand{\model}{\mathcal{M}} \]
Summary
In this problem set, we’re going to get practice using some new regression models to do the kinds of analysis we’ve been doing all semester. While we do this, we’ll pay a bit of attention to bias—a new phenomenon, because it’s not something that happens when we use the ‘all curves’ model we’ve been using.
Here are the libraries we’ll need.
And some settings to configure how our plots and code output looks.
theme_set(theme_minimal())
options(scipen = 1, digits = 3)Bone Density
Researchers hypothesize that higher calcium intake causes greater bone density (stronger bones). To test this, a clinic runs an experiment where patients are treated with a calcium supplement or a placebo for 1 year, after which their bone density is measured. Bone Density T-Scores indicate how good someone’s bone density is—0 is typical, positive numbers are better, negative numbers are worse (a threshold score of -2.5 is used to diagnose osteoporosis). The (fake) data from this imaginary study is shown below.
To get a random sample from its patient population in which treatment is randomized, the following procedure is used.
- Patients in our population are classified as old (age 65+) or young (age 64-). We’ll call this \(X\): \(X=0\) if young, \(X=1\) if old.
- For each patient in our population, we flip a coin to assign them to treatment or placebo. An age-group-specific weighted coin, so that patients are assigned treatment or placebo with probability depending on their age group. We’ll call this \(W\): \(W=1\) if treatment, \(W=0\) if placebo.
- We then sample with replacement from the population.1
Here’s a table of the data.
bone.data = data.frame(age = c(1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0),
supplement = c(1, 1, 0, 1, 1, 1, 0, 0, 0, 1, 0, 1),
bone_density = c(-2, -3, -2, 0, -1, 0, -2, 0, 1, 1, -1, 1.5))
bone.data| age | supplement | bone_density |
|---|---|---|
| 1 | 1 | -2.0 |
| 1 | 1 | -3.0 |
| 1 | 0 | -2.0 |
| 1 | 1 | 0.0 |
| 1 | 1 | -1.0 |
| 1 | 1 | 0.0 |
| 1 | 0 | -2.0 |
| 1 | 0 | 0.0 |
| 0 | 0 | 1.0 |
| 0 | 1 | 1.0 |
| 0 | 0 | -1.0 |
| 0 | 1 | 1.5 |
And here’s a plot. We’ll show the individual points, the mean and standard deviation of bone density for each age/treatment group (as points with arms), and the proportion of patients in each group that received treatment (as bar heights). Because we’re used to seeing data in which \(y\) is non-negative, we’ll plot bone density plus 3, so that our plot looks a bit more familiar.
plotX = function(w,x) { x + 0.15*(w-0.5) }
bone.plot =
ggplot(bone.data) +
geom_point(aes(x=plotX(supplement, age), y=bone_density+3, color=factor(supplement)),
alpha=1, size=2, position=position_jitter(width=0.05, height=0, seed=2)) +
stat_summary(aes(x=plotX(supplement, age), y=bone_density+3, color=factor(supplement)),
geom='pointrange', fun=mean,
fun.min=\(x)mean(x)-sd(x),
fun.max=\(x)mean(x)+sd(x), alpha=.5) +
geom_bar(aes(x=plotX(supplement, age), y=4*after_stat(prop), fill=factor(supplement), color=factor(supplement)),
position='identity', alpha=.1, linewidth=.2) +
scale_x_continuous(breaks=c(0,1), labels=c('Young', 'Old')) + xlab('')
bone.plot
Let’s start with some exercises where we estimate population summaries by taking appropriate averages of subsample means. These should be familiar from the last problem set, although there’s more of a causal flavor to their interpretation.
Now let’s recall an interpretation of what we’ve been doing Lecture 14. Using the subsample mean \(\hat\mu(w,x)\) to estimate the corresponding subpopulation mean \(\mu(w,x)\) is, equivalently, estimating the function \(\mu\) by the the function \(\hat\mu\) that minimizes the sum of squared residuals among all functions outright, i.e., \[ \hat \mu = \argmin_{m \in \model} \sum_{i=1}^n \qty{ Y_i - m(W_i,X_i) }^2 \qfor \model = \qty{ \text{ all functions } \ \ m(w,x) } \]
And we can, if we prefer, use the least squares estimator within a smaller set of functions, e.g., \[ \begin{aligned} \model &= \qty{ \text{ all functions of \emph{coarsened covariates} } \ \ m(w,x)=f(c(w,x)) } && \text{ is a piecewise-constant model } \\ \model &= \qty{ m(w,x) = m_0(w) + m_1(x) \ \ : \ \ m_0, m_1 \text{ are functions of one variable } } && \text{ is the additive model } \\ \model &= \qty{ m(w,x) = a(w) + b x \ \ : \ \ a \text{ is a function of one variable and } \ \ b \\ \text{ is a number } } && \text{ is the parallel lines model }\\ \model &= \qty{ m(w,x) = a(w) \ \ : \ \ a \text{ is a function of one variable } \ } && \text{ is the horizontal lines model }\\ \end{aligned} \]
Here’s code to find and plot you the least squares estimator \(\hat\mu\) within each model. To choose the coarsening \(c(w,x)\) in the piecewise-constant model, we’ll use a tree as implemented in the rpart package.
Here’s some generic code for making predictions and plotting given what R calls a ‘model fit object’.
prediction.function = function(model.fit) { function(w,x) { predict(model.fit, newdata=data.frame(supplement=w, age=x)) } }
plot.predictions = function(muhat) {
grid = expand.grid(supplement=c(0,1), age=c(0,1))
bone.plot +
geom_point(aes(x=age, y=muhat(supplement, age)+3, color=factor(supplement)), data=grid) +
geom_line(aes(x=age, y=muhat(supplement, age)+3, color=factor(supplement)), data=grid) +
guides(color='none', fill='none')
}And here’s some code to actually fit the models and plot the fits. The rpart stuff is its own thing, but for the rest of them, we’re using the R built-in function lm with an appropriate regression formula describing our model. The formula language takes a bit of getting used to and we’re not going to go get into that here. But if you’re curious, feel free to come to hours and ask about it.
## The piecewise-constant model
# Tell rpart to not to worry about having more than 1 observation per level of the coarsened covariate
# (or, in more evocative terms, per 'bucket')
rpart.settings = rpart.control(minbucket=1)
# Then fit and plot
piecewise.constant.model.fit = rpart(bone_density ~ supplement + age, data=bone.data, control=rpart.settings)
muhat.piecewise.constant = prediction.function(piecewise.constant.model.fit)
plot.predictions(muhat.piecewise.constant) + guides(color='none')
## The additive model
additive.model.fit = lm(bone_density ~ factor(supplement) + factor(age), data=bone.data)
muhat.additive = prediction.function(additive.model.fit)
plot.predictions(muhat.additive)
## The parallel lines model
parallel.lines.model.fit = lm(bone_density ~ factor(supplement) + age, data=bone.data)
muhat.parallel = prediction.function(parallel.lines.model.fit)
plot.predictions(muhat.parallel)
## The horizontal lines model
horizontal.lines.model.fit = lm(bone_density ~ factor(supplement), data=bone.data)
muhat.horizontal = prediction.function(horizontal.lines.model.fit)
plot.predictions(muhat.horizontal)
Note
If you look at these plots closely, you can see that the least squares estimators in the additive model and the parallel lines model are the same. In fact, when \(w\) and \(x\) are binary, these two models are the same—i.e., they’re two ways of describing the same set of functions. If you like, try to prove it.
Now let’s use a few of these models and take a moment to think about what happens.
Unless you actually believe these different summaries of the population are the same, this seems like a problem. And although our sample size is small, it’s not a small-sample-size problem. Using these models, we’ll get identical estimates of these different summaries no matter what sample size we have. This means at least one of the matching summary estimates has to be biased. This might be a sign that we’re using the wrong model and, in a sense, motivates our decision to use the ‘all functions’ model that we started with.
But before we move on, I want to call your attention to something. Unless you have a bug in your code (this is a good time to check), the raw difference in means you calculated in Exercise 3 and the corresponding part of Exercise 5 should be the same for every model. Let’s think about why that’s happening.
To make sense of this phenomenon, we’re going to think about the residuals we get when we predict \(Y_i\) using two functions in our model: the least squares estimator \(\hat \mu\) and a ‘perturbed version’ \(\hat\mu + t m\) for any other function \(m\) in our model. Or, to be precise, we’re going to think about the residuals we get when we predict \(Y_i\) using functions like this for a range of values of \(t\), i.e., for functions
\[ m_t(w,x) = \hat \mu(w,x) + t m(w,x) \quad \text{ for numbers } \ \ t \ \ \text{ near zero } \]
We know that the least squares estimator \(\hat \mu\) is the curve that minimizes the sum of squared residuals among all functions in our model, so we know that—if we think of the sum of squared residuals as a function of \(t\) for curves like this—the minimum occurs at \(t=0\) where \(m_t=\hat\mu\). Let’s see what the ‘zero-derivative’ condition characterizing this minimum tells us.
\[ \begin{aligned} 0 &= \frac{d}{dt} \sum_{i=1}^n \qty{ Y_i - m_t(W_i,X_i) }^2 \Big|_{t=0} && \text{ is the zero-derivative condition} \\ &= \sum_{i=1}^n 2\qty{ Y_i - m_\delta(W_i,X_i) } \times -\frac{d}{dt} m_{t}(W_i,X_i) \Big|_{t=0} && \text{ using the chain rule } \\ &= -2 \sum_{i=1}^n \qty{ Y_i - \hat\mu(W_i,X_i) } m(W_i,X_i) && \text{ differentiating } \ m_{t} \ \text{ and evaluating at } \ t=0 \\ \end{aligned} \]
This tells us that, for any function \(m\) in our model, if we multiply the least squares residuals \(\hat\varepsilon_i=Y_i-\hat\mu(W_i,X_i)\) by the corresponding values \(m(W_i,X_i)\) of that function and sum them up, we get zero. This is a generalization of the equation we used to show that the sample mean was the least squares estimator here in Lecture 14, where there was no multiplication by \(m(W_i,X_i)\)—where we just summed up the residuals.
Note
In linear algebra terms, this says that the vector \(\hat\varepsilon\) of least squares residuals is orthogonal to every vector of predictions \([m(W_1,X_1) \ldots m(W_n,X_n)]\) we can get using curves in our model. If you think of not orthogonal as meaning ‘looks a bit like’, this makes intuitive sense. If the residuals looked a bit like the predictions of some function \(m\) in our model, we’d be able to improve on the fit of \(\hat\mu\) by adding a bit of that function to it: \(\hat\mu + t m\) would, for small enough \(t\), fit the data better than \(\hat\mu\) would.
Now let’s think about what this tells us when we take \(m\) to be the constant function \(m(w,x)=1\) — assuming we’re in a model which, like every one we’ve discussed so far, includes that function. It tells us that the least squares estimator gets the overall mean right. \[ 0 = \sum_{i=1}^n \qty{ Y_i - \hat\mu(W_i,X_i) } \times 1 \quad \text{ so } \quad \frac{1}{n}\sum_{i=1}^n \hat\mu(W_i,X_i) = \frac{1}{n}\sum_{i=1}^n Y_i \]
That is, it tells us that the average of our predictions \(\hat\mu(W_i,X_i)\) over the people in our sample is the same as the average of their outcomes \(Y_i\). This is true for any model that includes the constant function \(m(w,x)=1\): the all-functions model, the additive model, the parallel lines model, the horizontal lines model, and most other models you can think of. It’s not true for the ‘lines through the origin’ model \(\{ m(w,x) = b(w)x \}\), which is why people really dislike that model.
We can say something similar about the average of our predictions \(\hat\mu(W_i,X_i)\) for people with a specific treatment status \(W_i=w\). If we use any model that includes an indicator for that treatment status, including any of the aforementioned four, then the average of our predictions \(\hat\mu(W_i,X_i)\) over the people in our sample with treatment status \(W_i=w\) is the same as the average of their outcomes \(Y_i\). \[ \begin{aligned} 0 &= \sum_{i=1}^n \qty{ Y_i - \hat\mu(W_i,X_i) } \times m(W_i,X_i) \qfor m(W_i,X_i) = \begin{cases} 1 & \text{ if } W_i=w \\ 0 & \text{ if } W_i \neq w \end{cases} \\ &= \sum_{i:W_i=w} \qty{ Y_i - \hat\mu(w,X_i) } \quad \text{ so } \frac{1}{N_w}\sum_{i:W_i=w} \hat\mu(w,X_i) = \frac{1}{N_w}\sum_{i:W_i=w} Y_i \qfor N_w=\sum_{i:W_i=w} 1 \end{aligned} \]
Tip
Recall that the difference in means we’re talking about looks like this.
\[ \frac{1}{N_w}\sum_{i:W_i=w} \hat\mu(w,X_i) - \frac{1}{N_{w'}}\sum_{i:W_i=w'} \hat\mu(w',X_i) \quad \text{for} \quad w=1 \ \text{ and } \ w'=0 \]
What do we know about each of the terms in this difference?
Gerber, Green, and Larimer Michigan Voter Experiment
Now let’s move on to a different dataset. This one is from a field experiment conducted by Gerber, Green, and Larimer in Michigan—the same one we discussed in Lecture 6.
Let’s start by loading the data. \(Y\) will be an indicator for whether the voter turned out in the 2006 election, \(X\) be the voter’s age at the time of the election, and \(W\) be an indicator for whether the voter received the treatment. \(W=1\) indicates that they received the treatment—a letter encouraging them to turn out—and \(W=0\) indicates that they did not. To make the data a bit easier to work with, I’m going to drop voters who are 96+ years old. What’s nice about this is that it ensures we’ve got voters in both treatment groups at every remaining age.
social.pressure = read.csv('social-pressure-data.csv')
keep = social.pressure$treatment %in% c(' Control', ' Neighbors') & (2006 - social.pressure$yob)< 96
turnout.data = data.frame(W=as.numeric(social.pressure[keep,'treatment'] == ' Neighbors'),
X=2006-social.pressure[keep, 'yob'],
Y = as.numeric(social.pressure[keep,'voted'] == 'Yes'))Here’s a plot a lot like the one we made for the bone density data. We’ll show the individual points, the mean and standard deviation of turnout for each age/treatment group (as points with arms), and the proportion of patients in each group that received treatment (as bar heights). Here the raw data is a bit harder to look at for a few reasons.
- There are way more points.
- There are way more columns.
- The outcome is binary, so all the points in each column stack on top of each other.
To make the plot more readable, we’ll make the points transparent and jitter them a bit—both horizontally and vertically — so they’re not quite on top of each other. It helps, but it’s still a bit hard to look at the raw data. The means and proportions are still easy enough to see, though. I’m also going to plot points for a random quarter of our survey recipients, rather than all of them, to help reduce the size of the images and therefore the strain they put on our web browser. This won’t have much of an effect visually.
plotX = function(w,x) { x + 0.75*(w-0.5) }
turnout.summary.plot =
ggplot(turnout.data) +
stat_summary(aes(x=plotX(W,X), y=Y, color=factor(W)),
geom='pointrange', fun=mean,
fun.min=\(x)mean(x)-sd(x),
fun.max=\(x)mean(x)+sd(x), alpha=.2, size=.5, linewidth=.5) +
geom_bar(aes(x=plotX(W, X), y=10*after_stat(prop), fill=factor(W), color=factor(W)),
position='identity', alpha=.1, linewidth=.1) +
guides(fill='none', color='none') + xlab('')
plot.point = sample(1:nrow(turnout.data), round(nrow(turnout.data)/4))
turnout.plot = turnout.summary.plot +
geom_point(aes(x=plotX(W, X), y=Y, color=factor(W)),
alpha=.01, size=.01, position=position_jitter(width=0.25, height=.1, seed=2), data=turnout.data[plot.point,])
turnout.plot 